1 Questions

Question 1: How is grazer density related to macrohabitat composition?

Question 2: How are vegetation structure and functional traits influenced by grazing?

2 Data preparation

#packages
library(sp)
library(sf)
## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
library(fields)
## Loading required package: spam
## Spam version 2.10-0 (2023-10-23) is loaded.
## Type 'help( Spam)' or 'demo( spam)' for a short introduction 
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix '.spam' to the function name, e.g. 'help( chol.spam)'.
## 
## Attaching package: 'spam'
## The following objects are masked from 'package:base':
## 
##     backsolve, forwardsolve
## Loading required package: viridisLite
## 
## Try help(fields) to get started.
library(lubridate)
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
library(ggplot2)
library(lattice)
library(reshape2)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(knitr)
library(arm)
## Loading required package: MASS
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## The following object is masked from 'package:spam':
## 
##     det
## Loading required package: lme4
## 
## arm (Version 1.14-4, built: 2024-4-1)
## Working directory is /Users/lillalovasz/Library/CloudStorage/Dropbox/papers/vegetation-grazer/rcode
## 
## Attaching package: 'arm'
## The following object is masked from 'package:spam':
## 
##     display
#read data
dhorse19_20 <- read.table("../data/horses-all-positions_2019-jan-09_2020_aug-31.txt", header=TRUE, sep = "\t")
dhorse20_21 <- read.table("../data/horses-all-positions_2020-sept-01_2021_aug_31.txt", header=TRUE, sep = "\t")
dcattle <- read.table("../data/cattle_all-positions-2019-jan-09_2021-aug-31.txt", header=TRUE, sep = "\t")
dfencefull <- read.csv("../data/fence.csv", header=TRUE)
dfence <- read.table("../data/fence-reduced.txt", header=TRUE, sep="\t")
dfence <- dfencefull # when we include the forest
dcells <- read.table("../data/cells_info_veg-plots-with-forest.txt", header=TRUE, sep="\t")
dlandcover <- read.table("../data/habitat-with-forest.txt",header=TRUE, sep="\t")

dveg <- read.table("../data/veg_data_2018-2019-2020-2021-summer.txt",header=TRUE, sep="\t")

2.1 prepare data for the grazer-habitat use analyses

Make all coordinates being the same
Clip grazer positions to the fence-area including the forest Correct the number of locations for the duration of the season and nr of animals tracked.

##horse data:
dhorse <- rbind(dhorse19_20,dhorse20_21)
head(dhorse)
##                                   Name     ID     Date  Time  X Latitude
## 1                              Onissoi 401776 1/9/2019 00:00 NA 47.62875
## 2                               Odissa 401777 1/9/2019 00:00 NA 47.62950
## 3                  TERMINATED - 401775 401775 1/9/2019 00:00 NA 47.62876
## 4                              Odaline 401778 1/9/2019 00:00 NA 47.62933
## 5 TERMINATED - 401774 lost (ex-Olocen) 401774 1/9/2019 00:00 NA 47.62883
## 6                              Odaline 401778 1/9/2019 01:00 NA 47.62976
##   Longitude Speed
## 1  7.556987  0.00
## 2  7.556098  0.00
## 3  7.556587  0.37
## 4  7.556327  0.00
## 5  7.556852  2.04
## 6  7.556088  0.74
str(dhorse)
## 'data.frame':    102216 obs. of  8 variables:
##  $ Name     : chr  "Onissoi" "Odissa" "TERMINATED - 401775" "Odaline" ...
##  $ ID       : int  401776 401777 401775 401778 401774 401778 401776 401774 401777 401775 ...
##  $ Date     : chr  "1/9/2019" "1/9/2019" "1/9/2019" "1/9/2019" ...
##  $ Time     : chr  "00:00" "00:00" "00:00" "00:00" ...
##  $ X        : logi  NA NA NA NA NA NA ...
##  $ Latitude : num  47.6 47.6 47.6 47.6 47.6 ...
##  $ Longitude: num  7.56 7.56 7.56 7.56 7.56 ...
##  $ Speed    : num  0 0 0.37 0 2.04 0.74 0.19 2.22 0.56 0.19 ...
dhorse$ID <- as.factor(dhorse$ID) 
dhorse$Name <- as.factor(dhorse$Name) 
levels(dhorse$Name) 
## [1] ""                                    
## [2] "Odaline"                             
## [3] "Odissa"                              
## [4] "Onissoi"                             
## [5] "TERMINATED - 401707"                 
## [6] "TERMINATED - 401774 lost (ex-Olocen)"
## [7] "TERMINATED - 401775"
dhorse$X <- NULL

length(which(is.na(dhorse)))
## [1] 27244
dhorse <- dhorse[complete.cases(dhorse), ] 

##cattle data:
head(dcattle,2)
##                   Name     ID     Date  Time  X Latitude Longitude X.1 Speed
## 1  TERMINATED - 401773 401773 1/9/2019 06:00 NA 47.62071  7.539350  NA  0.74
## 2  TERMINATED - 401771 401771 1/9/2019 06:01 NA 47.62059  7.538698  NA  0.19
dcattle$ID <- as.factor(dcattle$ID)
dcattle$Name <- as.factor(dcattle$Name)
levels(dcattle$ID)
## [1] "401767"   "401769"   "401770"   "401771"   "401773"   "421430"   "421431"  
## [8] "T5H-6839" "T5H-6840"
levels(dcattle$Name)
## [1] " TERMINATED - 401773"   "Gevaudan (6759) Tellus" "Kayla (3412) Tellus"   
## [4] "TERMINATED - 401767"    "TERMINATED - 401769"    "TERMINATED - 401770"   
## [7] "TERMINATED - 401771"    "TERMINATED - 421430"    "TERMINATED - 421431"
dcattle$X <- NULL
dcattle$X.1 <- NULL

##
#merge horse and cattle
dgrazer <- rbind(dhorse, dcattle)

#create "Animal" type 
dgrazer <- mutate(dgrazer, Animal = ifelse(dgrazer$ID %in% c("401707","401774","401775","401776","401777","401778"), "Horse", "Cattle"))
head(dgrazer,2)
##      Name     ID     Date  Time Latitude Longitude Speed Animal
## 1 Onissoi 401776 1/9/2019 00:00 47.62875  7.556987     0  Horse
## 2  Odissa 401777 1/9/2019 00:00 47.62950  7.556098     0  Horse
length(which(is.na(dgrazer)))
## [1] 0
dgrazer$Animal <- as.factor(dgrazer$Animal)

#date in correct format
index <- grepl("/", dgrazer$Date)
dgrazer$Date[index] <- as.character(as.Date(as.character(dgrazer$Date[index]), format="%m/%d/%Y"))
dgrazer$Date[!index] <- as.character(as.Date(as.character(dgrazer$Date[!index]), format="%d.%m.%y"))

#date for creating yday (day of year as decimal nr)
dgrazer$DateTime <- as.POSIXct(paste(as.character(dgrazer$Date), as.character(dgrazer$Time), sep=" "),format="%Y-%m-%d %H:%M", tz="CET")

##linear timeline as a number passed from 2019-jan-1.
dgrazer$dayofyear <- strptime(dgrazer$Date, format="%Y-%m-%d")$yday
dgrazer$year <- strptime(dgrazer$Date, format="%Y-%m-%d")$year+1900 # 
head(dgrazer,2)
##      Name     ID       Date  Time Latitude Longitude Speed Animal   DateTime
## 1 Onissoi 401776 2019-01-09 00:00 47.62875  7.556987     0  Horse 2019-01-09
## 2  Odissa 401777 2019-01-09 00:00 47.62950  7.556098     0  Horse 2019-01-09
##   dayofyear year
## 1         8 2019
## 2         8 2019
length(unique(dgrazer$dayofyear)) 
## [1] 366
#coordinates in correct format
tmp <- dgrazer
coordinates(tmp) <- c("Longitude", "Latitude") # make into a spatial point dataframe (WGS84 = long /y/; lat /x/)
crs.geo <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
proj4string(tmp) <- crs.geo #define projection
tmp2 <- spTransform(tmp, CRS("+init=epsg:32232")) # utm transformation
## Warning in CPL_crs_from_input(x): GDAL Message 1: +init=epsg:XXXX syntax is
## deprecated. It might return a CRS with a non-EPSG compliant axis order.
#dgrazer$Latitude.UTM <- coordinates(tmp2)[,"Latitude"]
#dgrazer$Longitude.UTM <- coordinates(tmp2)[,"Longitude"]
#
dgrazer$Latitude.UTM <- coordinates(tmp2)[,"coords.x2"]
dgrazer$Longitude.UTM <- coordinates(tmp2)[,"coords.x1"]


tmp3 <- dfence
coordinates(tmp3) <- c("x", "y") # make into a spatial point dataframe 
crs.geo <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
proj4string(tmp3) <- crs.geo #define projection
tmp4 <- spTransform(tmp3, CRS("+init=epsg:32232")) # utm transformation
dfence$Latitude.UTM <- coordinates(tmp4)[,"coords.x2"]
dfence$Longitude.UTM <- coordinates(tmp4)[,"coords.x1"]

plot(dfence$x, dfence$y)
text(dfence$x, dfence$y, labels=dfence$Name, cex=0.6)

# 
plot(dfence$Longitude.UTM, dfence$Latitude.UTM)

# text(dfence$Longitude.UTM, dfence$Latitude.UTM, labels=dfence$Name, cex=0.6)
head(dfence)
##       Name        y        x Latitude.UTM Longitude.UTM
## 1 Marker 1 47.62670 7.559454      5275811      391758.2
## 2 Marker 2 47.62683 7.559599      5275826      391769.4
## 3 Marker 3 47.62706 7.559288      5275852      391746.5
## 4 Marker 4 47.62732 7.558995      5275881      391725.1
## 5 Marker 5 47.62760 7.558683      5275913      391702.2
## 6 Marker 6 47.62788 7.558317      5275945      391675.3
########
#track and delete duplications
dgrazer$DateTime <- paste(as.character(dgrazer$Date), dgrazer$Time, dgrazer$Name)
table(table(as.character(dgrazer$DateTime)))
## 
##      1      2      3      4      5 
## 147695    363     79     11      5
index2 <- duplicated(dgrazer$DateTime)
dgrazer <- dgrazer[!index2,]

#plot raw grazer points (outside the grid too)
plot(dcells$x, dcells$y, type="n", xlim=c(min(dcells$x), max(dcells$x)+50),
     ylim=c(min(dcells$y), max(dcells$y)+50))
rect(dcells$x, dcells$y, dcells$x+50, dcells$y+50, border = "grey")
text(dcells$x, dcells$y, labels=dcells$cellnr, cex=0.6)

points(dgrazer$Longitude.UTM, dgrazer$Latitude.UTM, col=alpha(ifelse(dgrazer$Animal=='Horse', "blue","red"), 0.3), cex=0.6, pch=16)
points(dfence$Longitude.UTM, dfence$Latitude.UTM)

#points(dgrazer$Longitude.UTM, dgrazer$Latitude.UTM, col="orange", cex=0.6, pch=16)
 
#only positions inside the fence

dgrazer$Inside <- point.in.polygon(dgrazer$Longitude.UTM, dgrazer$Latitude.UTM, dfence$Longitude.UTM, dfence$Latitude.UTM)

dgrazer <- dgrazer[dgrazer$Inside==1,]

#all-fence-area
#plot animal position in all the enclosure 

plot(dfence$Longitude.UTM, dfence$Latitude.UTM)
points(dgrazer$Longitude.UTM, dgrazer$Latitude.UTM, col=alpha(ifelse(dgrazer$Animal=='Horse', "blue","red"), 0.3), cex=0.6, pch=16)

# grazer positions inside the grid cells 

allexistingx <- unique(dcells$x) 
allexistingy <- unique(dcells$y) 

for(i in 1:nrow(dgrazer)){
  index <- (dgrazer$Longitude.UTM[i]-allexistingx) >0 &  (dgrazer$Longitude.UTM[i]-allexistingx) <=50
  if(sum(index)>0) dgrazer$x.cell[i] <- allexistingx[index]
  index <- (dgrazer$Latitude.UTM[i]-allexistingy) >0 &  (dgrazer$Latitude.UTM[i]-allexistingy) <=50
  if(sum(index)>0) dgrazer$y.cell[i] <- allexistingy[index]
} 

dgrazer <- dgrazer[!is.na(dgrazer$y.cell),] 
#get the cellnr in the dataframe
dgrazer$cellnr <- dcells$cellnr[match(paste(dgrazer$x.cell, dgrazer$y.cell), paste(dcells$x, dcells$y))] 
length(unique(dgrazer$cellnr)) 
## [1] 160
# 

dgrazer <- dgrazer[!is.na(dgrazer$cellnr),] # -> delete positions outside the grid

###replot
plot(dcells$x, dcells$y, type="n", xlim=c(min(dcells$x), max(dcells$x)+50),
     ylim=c(min(dcells$y), max(dcells$y)+50))
rect(dcells$x, dcells$y, dcells$x+50, dcells$y+50, border = "grey")
text(dcells$x, dcells$y, labels=dcells$cellnr, cex=0.6)

points(dgrazer$Longitude.UTM, dgrazer$Latitude.UTM, cex=0.6,
       pch=16, col=rgb(1,1,0,0.2))

Question 1 (macrohabitat use of grazers): we use the average number of locations within each cell per season and year, averaged over the individuals tracked and averaged per day

dgrazer$Date.date <- strptime(dgrazer$Date, format="%Y-%m-%d")
dgrazerforveg <- dgrazer # make a copy for the micro-analyses later
dgrazer$season[dgrazer$Date.date >= "2019-06-24" & dgrazer$Date.date <= "2019-08-31"] <- "summer"
dgrazer$year.season[dgrazer$Date.date >= "2019-06-24" & dgrazer$Date.date <= "2019-08-31"] <- 2019
dgrazer$duration[dgrazer$Date.date >= "2019-06-24" & dgrazer$Date.date <= "2019-08-31"] <- 7+31+31

dgrazer$season[dgrazer$Date.date >= "2020-06-01" & dgrazer$Date.date <= "2020-08-31"] <- "summer"
dgrazer$year.season[dgrazer$Date.date >= "2020-06-01" & dgrazer$Date.date <= "2020-08-31"] <- 2020
dgrazer$duration[dgrazer$Date.date >= "2020-06-01" & dgrazer$Date.date <= "2020-08-31"] <- 30+31+31

dgrazer$season[dgrazer$Date.date >= "2021-06-28" & dgrazer$Date.date <= "2021-08-31"] <- "summer"
dgrazer$duration[dgrazer$Date.date >= "2021-06-28" & dgrazer$Date.date <= "2021-08-31"] <- 3+31+31
dgrazer$year.season[dgrazer$Date.date >= "2021-06-28" & dgrazer$Date.date <= "2021-08-31"] <- 2021


dgrazer$season[dgrazer$Date.date >= "2019-01-09" & dgrazer$Date.date <= "2019-02-28"] <- "winter"
dgrazer$duration[dgrazer$Date.date >= "2019-01-09" & dgrazer$Date.date <= "2019-02-28"] <- 23+28
dgrazer$year.season[dgrazer$Date.date >= "2019-01-09" & dgrazer$Date.date <= "2019-02-28"] <- 2019

dgrazer$season[dgrazer$Date.date >= "2019-12-01" & dgrazer$Date.date <= "2020-02-29"] <- "winter"
dgrazer$duration[dgrazer$Date.date >= "2019-12-01" & dgrazer$Date.date <= "2020-02-29"] <- 
31+31+29
dgrazer$year.season[dgrazer$Date.date >= "2019-12-01" & dgrazer$Date.date <= "2020-02-29"] <- 2020

dgrazer$season[dgrazer$Date.date >= "2020-12-01" & dgrazer$Date.date <= "2021-02-28"] <- "winter"
dgrazer$duration[dgrazer$Date.date >= "2020-12-01" & dgrazer$Date.date <= "2021-02-28"] <- 31+31+28
dgrazer$year.season[dgrazer$Date.date >= "2020-12-01" & dgrazer$Date.date <= "2021-02-28"] <- 2021

dgrazer <- dgrazer[!is.na(dgrazer$season),]

dgrazer$year.season <- factor(dgrazer$year.season)
dgrazer$season <- factor(dgrazer$season)
dgrazer$cellnr <- factor(dgrazer$cellnr)

dgrazersum <- aggregate(Animal~year.season+season+cellnr, dgrazer, FUN=function(x) sum(x=="Cattle"), drop=FALSE)
names(dgrazersum)[names(dgrazersum)=="Animal"] <- "Cattle"
dgrazersum$Horse <- aggregate(Animal~year.season+season+cellnr, dgrazer, FUN=function(x) sum(x=="Horse"), drop=FALSE)$Animal
dgrazersum$Cattle[is.na(dgrazersum$Cattle)] <- 0
dgrazersum$Horse[is.na(dgrazersum$Horse)] <- 0
dgrazersum$duration <- dgrazer$duration[match(paste(dgrazersum$year.season, dgrazersum$season), paste(dgrazer$year.season, dgrazer$season))]
dgrazersum$grazer <- dgrazersum$Cattle + dgrazersum$Horse

dgrazersum$grazer <- dgrazersum$grazer/dgrazersum$duration
dgrazersum$Cattle <- dgrazersum$Cattle/dgrazersum$duration
dgrazersum$Horse <- dgrazersum$Horse/dgrazersum$duration


# number of individuals tracked
dtemp <- aggregate(ID~year+season, dgrazer[dgrazer$Animal=="Horse",], FUN=function(x) length(unique(x)), drop=FALSE)
dgrazersum$Horse_nrId <- dtemp$ID[match(paste(dgrazersum$year, dgrazersum$season),
                                        paste(dtemp$year, dtemp$season))]
dtemp <- aggregate(ID~year+season, dgrazer[dgrazer$Animal=="Cattle",], FUN=function(x) length(unique(x)), drop=FALSE)
dgrazersum$Cattle_nrId <- dtemp$ID[match(paste(dgrazersum$year, dgrazersum$season),
                                        paste(dtemp$year, dtemp$season))]


dgrazersum$horseUDcorr <- dgrazersum$Horse/dgrazersum$Horse_nrId
dgrazersum$cattleUDcorr <- dgrazersum$Cattle/dgrazersum$Cattle_nrId

dgrazersum$grazerUDcorr <- dgrazersum$horseUDcorr+dgrazersum$cattleUDcorr
hist(dgrazersum$grazerUDcorr)

range(dgrazersum$grazerUDcorr)
## [1] 0.000000 3.813043
dgrazersum$x <- dcells$x[match(dgrazersum$cellnr, dcells$cellnr)]
dgrazersum$y <- dcells$y[match(dgrazersum$cellnr, dcells$cellnr)]
dgrazersum$grazer.log <- log(dgrazersum$grazerUDcorr+min(dgrazersum$grazer[dgrazersum$grazer>0])/2)
dgrazersum$horse.log <- log(dgrazersum$horseUDcorr++min(dgrazersum$Horse[dgrazersum$Horse>0])/2)
dgrazersum$cattle.log <- log(dgrazersum$cattleUDcorr++min(dgrazersum$Cattle[dgrazersum$Cattle>0])/2)

maxloggrazer <- max(dgrazersum$grazer.log)
minloggrazer <- min(dgrazersum$grazer.log)
dgrazersum$colnum <- (dgrazersum$grazer.log-minloggrazer)/(maxloggrazer-minloggrazer)

par(mfrow=c(3,2), mar=c(1,1,2,1))
for(i in levels(dgrazersum$year.season)){
  index <- dgrazersum$season=="winter"&dgrazersum$year.season==i 
  plot(dgrazersum$x[index], dgrazersum$y[index], pch=15, col=grey(1-dgrazersum$colnum[index]), cex=3,    main=paste("Winter",i))
  index <- dgrazersum$season=="summer"&dgrazersum$year.season==i
  plot(dgrazersum$x[index], dgrazersum$y[index], pch=15, col=grey(1-dgrazersum$colnum[index]), cex=3, main=paste("Summer", i))
}

##put landcover in the dataframe
#head(dlandcover)
dat <- merge(dgrazersum, dlandcover, by="cellnr")

#head(dat,2)
cov.varnames <- c("Cover.tree", "Cover.sapling", "Cover.shrub", "Cover.meadow", "Surface.water", "Bareground") 
pairs(dat[, cov.varnames])
Pairs-plot of the cover variables

Figure 2.1: Pairs-plot of the cover variables

2.2 prepare data for analysis of the plant traits on a microhabitat scale

Question 2 (how do vegetation traits change as a function of grazer density): we use the grazer density as a predictor for the vegetation structure. Therefore, we need to measure the total grazing density summed over all individuals. Because not all individuals were tracked, and because both grazer species are known to stick together in group, we multiplied the average locations per individual by the number of animals that were actually grazing.

# need to add in dveg the grazing density 

dgrazerforveg$include <- FALSE
dgrazerforveg$include[dgrazerforveg$Date.date >= "2019-01-09" & dgrazerforveg$Date.date <= "2019-06-01"] <- TRUE
dgrazerforveg$year.season[dgrazerforveg$Date.date >= "2019-01-09" & dgrazerforveg$Date.date <= "2019-06-01"] <- 2019

dgrazerforveg$include[dgrazerforveg$Date.date >= "2020-01-01" & dgrazerforveg$Date.date <= "2020-06-01"] <- TRUE
dgrazerforveg$year.season[dgrazerforveg$Date.date >= "2020-01-01" & dgrazerforveg$Date.date <= "2020-06-01"] <- 2020

dgrazerforveg$include[dgrazerforveg$Date.date >= "2021-01-01" & dgrazerforveg$Date.date <= "2021-06-01"] <- TRUE
dgrazerforveg$year.season[dgrazerforveg$Date.date >= "2021-01-01" & dgrazerforveg$Date.date <= "2021-06-01"] <- 2021

dgrazerforveg <- dgrazerforveg[dgrazerforveg$include,]

dgrazersumveg <- aggregate(Animal~year.season+cellnr, dgrazerforveg, FUN=function(x) sum(x=="Cattle"), drop=FALSE)
names(dgrazersumveg)[names(dgrazersumveg)=="Animal"] <- "Cattle"
dgrazersumveg$Horse <- aggregate(Animal~year.season+cellnr, dgrazerforveg, FUN=function(x) sum(x=="Horse"), drop=FALSE)$Animal
dgrazersumveg$Cattle[is.na(dgrazersumveg$Cattle)] <- 0
dgrazersumveg$Horse[is.na(dgrazersumveg$Horse)] <- 0

# number of individuals tracked measured in "animal-days"
dtemp<- aggregate(dayofyear~Name + year, data=dgrazerforveg, function(x) length(unique(x)))
dtemp$Animal <- dgrazer$Animal[match(dtemp$Name, dgrazer$Name)]
dtemp1 <- aggregate(dayofyear~year+Animal, data=dtemp, sum)
dgrazersumveg$nrCattletracked <- dtemp1$dayofyear[dtemp1$Animal=="Cattle"][match(dgrazersumveg$year.season, dtemp1$year[dtemp1$Animal=="Cattle"])]
dgrazersumveg$nrHorsetracked <- dtemp1$dayofyear[dtemp1$Animal=="Horse"][match(dgrazersumveg$year.season, dtemp1$year[dtemp1$Animal=="Horse"])]


drealpres <- read.table("../data/animal-real-presence.txt",  header=TRUE, sep="\t")
drealpres$start.date.date <- strptime(drealpres$start.date, format="%d/%m/%Y")
drealpres$end.date.date <- strptime(drealpres$end.date, format="%d/%m/%Y")
dtemp2 <- data.frame(date=seq(min(drealpres$start.date.date),max(drealpres$end.date.date), by=60*60*24))
dtemp2$dayofyear <- yday(dtemp2$date)
dtemp2 <- dtemp2[dtemp2$dayofyear<=yday(strptime("01.06.2020", format="%d.%m.%Y")),]

for(i in 1:sum(drealpres$animal=="horse")){
  ii <- c(1:nrow(drealpres))[drealpres$animal=="horse"][i]
  dtemp2$horse[dtemp2$date>=drealpres$start.date.date[ii] & dtemp2$date<=drealpres$end.date.date[ii]] <- drealpres$number[ii]
}
for(i in 1:sum(drealpres$animal=="cattle")){
  ii <- c(1:nrow(drealpres))[drealpres$animal=="cattle"][i]
  dtemp2$cattle[dtemp2$date>=drealpres$start.date.date[ii] & dtemp2$date<=drealpres$end.date.date[ii]] <- drealpres$number[ii]
}
dtemp2$cattle[is.na(dtemp2$cattle)] <- 0
dtemp2$horse[is.na(dtemp2$horse)] <- 0
dtemp2$year <- year(dtemp2$date)
dtemp3 <- aggregate(cattle~year, data=dtemp2, FUN=sum)
dtemp3$horse <- aggregate(horse~year, data=dtemp2, FUN=sum)$horse

dgrazersumveg$realpresCattle <- dtemp3$cattle[match(dgrazersumveg$year.season, dtemp3$year)]
dgrazersumveg$realpresHorse <- dtemp3$horse[match(dgrazersumveg$year.season, dtemp3$year)]
dgrazersumveg$cattle.corr <- dgrazersumveg$Cattle/dgrazersumveg$nrCattletracked*dgrazersumveg$realpresCattle
dgrazersumveg$horse.corr <-dgrazersumveg$Horse/dgrazersumveg$nrHorsetracked*dgrazersumveg$realpresHorse
dgrazersumveg$grazer.corr <- dgrazersumveg$horse.corr + dgrazersumveg$cattle.corr

Information about the dataset Flora Indicativa (Landolt et al 2010), which is used to assign functional traits to the sampled plant species.

The following variables of the Flora Indicativa are used:

  • Climate - L (light): 1 deep shade; 2 shade; 3 semi-shade; 4 well lit; 5 full light [additionally we -> sd light because grazers could affect the between-patch variance]

  • Soil - N (nutrients): 1 very infertile; 2 infertile; 3 medium infertile to medium fertile; 4 fertile; 5 very fertile and over-rich

  • Life history - MV (mowing /grazing tolerance) 1 bad; 2 little; 3 moderate; 4 good; 5 very good

(the variable name is in two rows in the original excel file)

Information about the dataset of field plant sampling: The raw dataset contains typos, which will be corrected below (the original txt dataset remains unchanged). Elymus pungens occurs in the Petite Camargue but no data exist in Flora Helvetica.
Salix x rubens is a hybrid species in the genus Salix. As the hybrid is not present in the Flora Indicativa, we use the attributes of the type species of the genus, Salix alba. Species that were not determined to the species level -> NA

## cellnr corresponds to plot_nr
dcellnr <- read.table("../data/cells_info_veg-plots-with-forest.txt", header=TRUE, sep="\t")
dveg$cellnr <- dcellnr$cellnr[match(dveg$plot_nr, dcellnr$plot_nr)]

dflora <- read.table("../data/FLORA-INDICATIVA-EN.txt", header=TRUE, skip=1, sep="\t")
dflora$L..light. <- as.numeric(dflora$L..light.) # x turn to NA
## Warning: NAs introduced by coercion
dveg$species1 <- trimws(dveg$species1, which="both")
dveg$species2 <- trimws(dveg$species2, which="both")
dveg$species3 <- trimws(dveg$species3, which="both")

# correct species names
index <- dveg$species1=="Calamagrostis epigeios"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Calamagrostis epigejos"
index <- dveg$species2=="Calamagrostis epigeios"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Calamagrostis epigejos"
index <- dveg$species3=="Calamagrostis epigeios"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Calamagrostis epigejos"


index <- dveg$species1=="Calamgrostis epigeios"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Calamagrostis epigejos"
index <- dveg$species2=="Calamgrostis epigeios"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Calamagrostis epigejos"
index <- dveg$species3=="Calamgrostis epigeios"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Calamagrostis epigejos"

index <- dveg$species1=="Seculigera varia"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Securigera varia"
index <- dveg$species2=="Seculigera varia"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Securigera varia"
index <- dveg$species3=="Seculigera varia"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Securigera varia"

index <- dveg$species1=="Seculigera Varia"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Securigera varia"
index <- dveg$species2=="Seculigera Varia"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Securigera varia"
index <- dveg$species3=="Seculigera Varia"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Securigera varia"

index <- dveg$species1=="Salix rubens"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Salix alba"
index <- dveg$species2=="Salix rubens"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Salix alba"
index <- dveg$species3=="Salix rubens"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Salix alba"

index <- dveg$species1=="Odontites serotinus"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Odontites vernus"
index <- dveg$species2=="Odontites serotinus"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Odontites vernus"
index <- dveg$species3=="Odontites serotinus"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Odontites vernus"

index <- dveg$species1=="Odontites rubra"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Odontites vernus"
index <- dveg$species2=="Odontites rubra"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Odontites vernus"
index <- dveg$species3=="Odontites rubra"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Odontites vernus"

index <- is.element(dveg$species1, 
                    c("Bromus hordaceous", "Bromus hordaeceus", "Bromus hordaceus", "Bromus hordaceous", "Bromes hordaceous")); index[is.na(index)] <- FALSE
dveg$species1[index] <- "Bromus hordeaceus" 
index <- is.element(dveg$species2, 
                    c("Bromus hordaceous", "Bromus hordaeceus", "Bromus hordaceus", "Bromus hordaceous", "Bromes hordaceous")); index[is.na(index)] <- FALSE
dveg$species2[index] <- "Bromus hordeaceus" 
index <- is.element(dveg$species3, 
                    c("Bromus hordaceous", "Bromus hordaeceus", "Bromus hordaceus", "Bromus hordaceous", "Bromes hordaceous")); index[is.na(index)] <- FALSE
dveg$species3[index] <- "Bromus hordeaceus" 

index <- dveg$species1=="Bromes sterilis"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Bromus sterilis"
index <- dveg$species2=="Bromes sterilis"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Bromus sterilis"
index <- dveg$species3=="Bromes sterilis"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Bromus sterilis"


index <- is.element(dveg$species1, 
                    c("Arenaria serpillyfolia", "Arenaria serpillifolia")); index[is.na(index)] <- FALSE
dveg$species1[index] <- "Arenaria serpyllifolia" 
index <- is.element(dveg$species2, 
                    c("Arenaria serpillyfolia", "Arenaria serpillifolia")); index[is.na(index)] <- FALSE
dveg$species2[index] <- "Arenaria serpyllifolia" 

index <- is.element(dveg$species3, 
                    c("Arenaria serpillyfolia", "Arenaria serpillifolia")); index[is.na(index)] <- FALSE
dveg$species3[index] <- "Arenaria serpyllifolia" 


index <- is.element(dveg$species1, 
                    c("Arrhenaterum elatius", "Arrhenatrum elatius")); index[is.na(index)] <- FALSE
dveg$species1[index] <- "Arrhenatherum elatius"
index <- is.element(dveg$species2, 
                    c("Arrhenaterum elatius", "Arrhenatrum elatius")); index[is.na(index)] <- FALSE
dveg$species2[index] <- "Arrhenatherum elatius"
index <- is.element(dveg$species3, 
                    c("Arrhenaterum elatius", "Arrhenatrum elatius")); index[is.na(index)] <- FALSE
dveg$species3[index] <- "Arrhenatherum elatius"

index <- dveg$species1=="Medicago campestre"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Trifolium campestre"
index <- dveg$species2=="Medicago campestre"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Trifolium campestre"
index <- dveg$species3=="Medicago campestre"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Trifolium campestre"

index <- dveg$species1=="Festuca  pratensis"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Festuca pratensis"
index <- dveg$species2=="Festuca  pratensis"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Festuca pratensis"
index <- dveg$species3=="Festuca  pratensis"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Festuca pratensis"

index <- is.element(dveg$species1, 
                    c("APetraraghia prolifera", "Petraraghia prolifera", "Petrarhaghia prolifera", "Petroraghia prolifera")); index[is.na(index)] <- FALSE
dveg$species1[index] <- "Petrorhagia prolifera"
index <- is.element(dveg$species2, 
                    c("APetraraghia prolifera", "Petraraghia prolifera", "Petrarhaghia prolifera", "Petroraghia prolifera")); index[is.na(index)] <- FALSE
dveg$species2[index] <- "Petrorhagia prolifera"
index <- is.element(dveg$species3, 
                    c("APetraraghia prolifera", "Petraraghia prolifera", "Petrarhaghia prolifera", "Petroraghia prolifera")); index[is.na(index)] <- FALSE
dveg$species3[index] <- "Petrorhagia prolifera"


index <- dveg$species1=="Gallium mollugo"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Galium mollugo"
index <- dveg$species2=="Gallium mollugo"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Galium mollugo"
index <- dveg$species3=="Gallium mollugo"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Galium mollugo"

index <- dveg$species1=="Plantage lanceolata"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Plantago lanceolata"
index <- dveg$species2=="Plantage lanceolata"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Plantago lanceolata"
index <- dveg$species3=="Plantage lanceolata"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Plantago lanceolata"

index <- dveg$species1=="Erigeron annus"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Erigeron annuus"
index <- dveg$species2=="Erigeron annus"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Erigeron annuus"
index <- dveg$species3=="Erigeron annus"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Erigeron annuus"

index <- dveg$species1=="Carex sp"; index[is.na(index)] <- FALSE
dveg$species1[index] <- NA
index <- dveg$species2=="Carex sp"; index[is.na(index)] <- FALSE
dveg$species2[index] <- NA
index <- dveg$species3=="Carex sp"; index[is.na(index)] <- FALSE
dveg$species3[index] <- NA


index <- is.element(dveg$species1, 
                    c("Juncus inflxus", "Juncos inflexus")); index[is.na(index)] <- FALSE
dveg$species1[index] <- "Juncus inflexus"
index <- is.element(dveg$species2, 
                    c("Juncus inflxus", "Juncos inflexus")); index[is.na(index)] <- FALSE
dveg$species2[index] <- "Juncus inflexus"
index <- is.element(dveg$species3, 
                    c("Juncus inflxus", "Juncos inflexus")); index[is.na(index)] <- FALSE
dveg$species3[index] <- "Juncus inflexus"

index <- dveg$species1=="Myosotis arvense"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Myosotis arvensis"
index <- dveg$species2=="Myosotis arvense"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Myosotis arvensis"
index <- dveg$species3=="Myosotis arvense"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Myosotis arvensis"

index <- dveg$species1=="Tussilagi farfara"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Tussilago farfara"
index <- dveg$species2=="Tussilagi farfara"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Tussilago farfara"
index <- dveg$species3=="Tussilagi farfara"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Tussilago farfara"

index <- dveg$species1=="Phragmites communis"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Phragmites australis"
index <- dveg$species2=="Phragmites communis"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Phragmites australis"
index <- dveg$species3=="Phragmites communis"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Phragmites australis"

index <- dveg$species1=="Argrostis stolonifera"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Agrostis stolonifera"
index <- dveg$species2=="Argrostis stolonifera"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Agrostis stolonifera"
index <- dveg$species3=="Argrostis stolonifera"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Agrostis stolonifera"

index <- dveg$species1=="Amorpha sp"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Amorpha fruticosa"
index <- dveg$species2=="Amorpha sp"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Amorpha fruticosa"
index <- dveg$species3=="Amorpha sp"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Amorpha fruticosa"

index <- dveg$species1=="Cirsium sp"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Cirsium vulgare"
index <- dveg$species2=="Cirsium sp"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Cirsium vulgare"
index <- dveg$species3=="Cirsium sp"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Cirsium vulgare"
# (this is the most common...)

index <- dveg$species1=="Echyum vulgare"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Echium vulgare"
index <- dveg$species2=="Echyum vulgare"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Echium vulgare"
index <- dveg$species3=="Echyum vulgare"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Echium vulgare"

index <- dveg$species1=="Menta aquatica"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Mentha aquatica"
index <- dveg$species2=="Menta aquatica"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Mentha aquatica"
index <- dveg$species3=="Menta aquatica"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Mentha aquatica"

index <- dveg$species1=="Medication lupulina"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Medicago lupulina"
index <- dveg$species2=="Medication lupulina"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Medicago lupulina"
index <- dveg$species3=="Medication lupulina"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Medicago lupulina"

index <- dveg$species1=="Cinosorus cristatus"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Cynosurus cristatus"
index <- dveg$species2=="Cinosorus cristatus"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Cynosurus cristatus"
index <- dveg$species3=="Cinosorus cristatus"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Cynosurus cristatus"

index <- dveg$species1=="Argostis gigantea"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Agrostis gigantea"
index <- dveg$species2=="Argostis gigantea"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Agrostis gigantea"
index <- dveg$species3=="Argostis gigantea"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Agrostis gigantea"

index <- dveg$species1=="Dactilys glomerata"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Dactylis glomerata"
index <- dveg$species2=="Dactilys glomerata"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Dactylis glomerata"
index <- dveg$species3=="Dactilys glomerata"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Dactylis glomerata"


index <- is.element(dveg$species1, 
                    c("Picris hieraceoides", "Picris hieraceous", "Picris hyeracioides")); index[is.na(index)] <- FALSE
dveg$species1[index] <- "Picris hieracioides"
index <- is.element(dveg$species2, 
                    c("Picris hieraceoides", "Picris hieraceous", "Picris hyeracioides")); index[is.na(index)] <- FALSE
dveg$species2[index] <- "Picris hieracioides"
index <- is.element(dveg$species3, 
                    c("Picris hieraceoides", "Picris hieraceous", "Picris hyeracioides")); index[is.na(index)] <- FALSE
dveg$species3[index] <- "Picris hieracioides"

index <- is.element(dveg$species1, 
                    c("Hypericum officinalis")); index[is.na(index)] <- FALSE
dveg$species1[index] <- "Hypericum perforatum"
index <- is.element(dveg$species2, 
                    c("Hypericum officinalis")); index[is.na(index)] <- FALSE
dveg$species2[index] <- "Hypericum perforatum"
index <- is.element(dveg$species3, 
                    c("Hypericum officinalis")); index[is.na(index)] <- FALSE
dveg$species3[index] <- "Hypericum perforatum"


index <- dveg$species1=="Potentilla reptans ?"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Potentilla reptans" 
index <- dveg$species2=="Potentilla reptans ?"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Potentilla reptans" 
index <- dveg$species3=="Potentilla reptans ?"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Potentilla reptans" 

#dveg[dveg==“Populus balsamifera”] <- NA 
# (Flora Indicativa does not contain, and not common in the PCA and in the data)

index <- dveg$species1=="Vicia cracca tenuifolia"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Vicia cracca"
index <- dveg$species2=="Vicia cracca tenuifolia"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Vicia cracca"
index <- dveg$species3=="Vicia cracca tenuifolia"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Vicia cracca"

index <- dveg$species1=="Origanium vulgare"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Origanum vulgare"
index <- dveg$species2=="Origanium vulgare"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Origanum vulgare"
index <- dveg$species3=="Origanium vulgare"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Origanum vulgare"

# end of correction of species names
dveg$Light_sp1 <- dflora$L..light.[match(dveg$species1, dflora$Taxon)]
dveg$Light_sp2 <- dflora$L..light.[match(dveg$species2, dflora$Taxon)]
dveg$Light_sp3 <- dflora$L..light.[match(dveg$species3, dflora$Taxon)]

unique(dveg$species1[is.na(dveg$Light_sp1)])
## [1] NA               "Elymus pungens" ""
unique(dveg$species2[is.na(dveg$Light_sp2)])
## [1] NA               "Elymus pungens" "Ranunculus sp"  "Artemisia sp"  
## [5] ""
unique(dveg$species3[is.na(dveg$Light_sp3)])
## [1] NA                    "Populus balsamifera" ""
dveg$N_sp1 <- dflora$N..Nutrients.[match(dveg$species1, dflora$Taxon)]
dveg$N_sp2 <- dflora$N..Nutrients.[match(dveg$species2, dflora$Taxon)]
dveg$N_sp3 <- dflora$N..Nutrients.[match(dveg$species3, dflora$Taxon)]

dveg$MV_sp1 <- dflora$MV..mowing.tolerance.[match(dveg$species1, dflora$Taxon)]
dveg$MV_sp2 <- dflora$MV..mowing.tolerance.[match(dveg$species2, dflora$Taxon)]
dveg$MV_sp3 <- dflora$MV..mowing.tolerance.[match(dveg$species3, dflora$Taxon)]


dveg$light.mean <- apply(dveg[,c("Light_sp1","Light_sp2","Light_sp3")], 1, mean, na.rm=TRUE)
dveg$light.sd <- apply(dveg[,c("Light_sp1","Light_sp2","Light_sp3")], 1, sd, na.rm=TRUE)

dveg$N_sp1 <- as.numeric(dveg$N_sp1)
dveg$N_sp2 <- as.numeric(dveg$N_sp2)
dveg$N_sp3 <- as.numeric(dveg$N_sp3)

dveg$N.mean <- apply(dveg[,c("N_sp1","N_sp2","N_sp3")], 1, mean, na.rm=TRUE)


dveg$MV_sp1 <- as.numeric(dveg$MV_sp1)
## Warning: NAs introduced by coercion
dveg$MV_sp2 <- as.numeric(dveg$MV_sp2)
## Warning: NAs introduced by coercion
dveg$MV_sp3 <- as.numeric(dveg$MV_sp3)
## Warning: NAs introduced by coercion
dveg$MV.mean <- apply(dveg[,c("MV_sp1","MV_sp2","MV_sp3")], 1, mean, na.rm=TRUE)


# add grazerdensity

dveg$horse.corr <- dgrazersumveg$horse.corr[match(paste(dveg$cellnr, dveg$year), paste(dgrazersumveg$cellnr, dgrazersumveg$year.season))]
dveg$cattle.corr <- dgrazersumveg$cattle.corr[match(paste(dveg$cellnr, dveg$year), paste(dgrazersumveg$cellnr, dgrazersumveg$year.season))]
dveg$grazer.corr <- dgrazersumveg$grazer.corr[match(paste(dveg$cellnr, dveg$year), paste(dgrazersumveg$cellnr, dgrazersumveg$year.season))]

dveg$grazer.corr[dveg$year==2018] <- 0
dveg$horse.corr[dveg$year==2018] <- 0
dveg$cattle.corr[dveg$year==2018] <- 0
dveg <- dveg[!is.na(dveg$year),]


dveg$heightm <- apply(dveg[,c("height1_1", "height1_2", "height1_3" ,"height1_4" ,  "height1_5" , "height1_6",  "height1_7" , "height1_8" ,  "height1_9"  ,
"height1_10" , "height2_1","height2_2",  "height2_3", "height2_4", "height2_5",
"height2_6", "height2_7", "height2_8", "height2_9", "height2_10")], 1, mean)


pairs(dveg[,c("grazer.corr", "horse.corr", "cattle.corr", "light.mean", "light.sd", "N.mean", "MV.mean", "heightm", "total_cover")], panel=panel.smooth)

3 Analysis

3.1 Grazer density in relation to cover variables on the macrohabitat scale

mod <- lmer(grazer.log~Cover.tree + Cover.sapling + Cover.shrub+ Surface.water+ Bareground + I(Cover.tree^2) + I(Cover.sapling^2) + I(Cover.shrub^2)+ I(Surface.water^2)+ I(Bareground^2) + season + season:(Cover.tree + Cover.sapling + Cover.shrub+ Surface.water+ Bareground + I(Cover.tree^2) + I(Cover.sapling^2) + I(Cover.shrub^2)+ I(Surface.water^2)+ I(Bareground^2)) + (1|cellnr), data=dat)
## Warning: Some predictor variables are on very different scales: consider
## rescaling
resdat <- data.frame(resid=resid(mod), x=dcells$x[match(dat$cellnr, dcells$cellnr)],
                     y=dcells$y[match(dat$cellnr, dcells$cellnr)])
variomod <- vgram(resdat[,c("x", "y")], resdat$resid, 
                    type="variogram")


par(mfrow=c(1,2))
plot(variomod, ylim=c(0,0.8))

# bubble plot
plot(jitter(dat$x,20), jitter(dat$y, 20), pch=16, col=c("blue", "orange")[(sign(resid(mod))+3)/2], cex=abs(resid(mod)))
Semi-variogram and bubble plot of the residuals of the grazer model to check for spatial correlation. Positive residuals are orange, negative residuals are blue.

Figure 3.1: Semi-variogram and bubble plot of the residuals of the grazer model to check for spatial correlation. Positive residuals are orange, negative residuals are blue. 

#plot(mod)
modh <- lmer(horse.log~Cover.tree + Cover.sapling + Cover.shrub+ Surface.water+ Bareground + I(Cover.tree^2) + I(Cover.sapling^2) + I(Cover.shrub^2)+ I(Surface.water^2)+ I(Bareground^2) + season + season:(Cover.tree + Cover.sapling + Cover.shrub+ Surface.water+ Bareground + I(Cover.tree^2) + I(Cover.sapling^2) + I(Cover.shrub^2)+ I(Surface.water^2)+ I(Bareground^2)) + (1|cellnr), data=dat)
## Warning: Some predictor variables are on very different scales: consider
## rescaling
resdat <- data.frame(resid=resid(modh), x=dcells$x[match(dat$cellnr, dcells$cellnr)],
                     y=dcells$y[match(dat$cellnr, dcells$cellnr)])
variomod <- vgram(resdat[,c("x", "y")], resdat$resid, 
                    type="variogram")


par(mfrow=c(1,2))
plot(variomod, ylim=c(0,0.8))

# bubble plot
plot(jitter(dat$x,20), jitter(dat$y, 20), pch=16, col=c("blue", "orange")[(sign(resid(modh))+3)/2], cex=abs(resid(mod)))
Visualisation of the residuals of the horse model to check for spatial correlation. Positive residuals are orange, negative residuals are blue.

Figure 3.2: Visualisation of the residuals of the horse model to check for spatial correlation. Positive residuals are orange, negative residuals are blue. 

#plot(modh)
modc <- lmer(cattle.log~Cover.tree + Cover.sapling + Cover.shrub+ Surface.water+ Bareground + I(Cover.tree^2) + I(Cover.sapling^2) + I(Cover.shrub^2)+ I(Surface.water^2)+ I(Bareground^2) + season + season:(Cover.tree + Cover.sapling + Cover.shrub+ Surface.water+ Bareground + I(Cover.tree^2) + I(Cover.sapling^2) + I(Cover.shrub^2)+ I(Surface.water^2)+ I(Bareground^2)) + (1|cellnr), data=dat)
## Warning: Some predictor variables are on very different scales: consider
## rescaling
resdat <- data.frame(resid=resid(modc), x=dcells$x[match(dat$cellnr, dcells$cellnr)],
                     y=dcells$y[match(dat$cellnr, dcells$cellnr)])
variomod <- vgram(resdat[,c("x", "y")], resdat$resid, 
                    type="variogram")


par(mfrow=c(1,2))
plot(variomod)

# bubble plot
plot(jitter(dat$x,20), jitter(dat$y, 20), pch=16, col=c("blue", "orange")[(sign(resid(modc))+3)/2], cex=abs(resid(mod)))
Visualisation of the residuals of the cattle model to check for spatial correlation. Positive residuals are orange, negative residuals are blue.

Figure 3.3: Visualisation of the residuals of the cattle model to check for spatial correlation. Positive residuals are orange, negative residuals are blue. 

#plot(modc)

No obvious spatial correlation can be seen in the residuals, i.e., both negative and positive residuals are spread around the whole study area. However, in the eastern part of the study area, the absolute values of the residuals are larger. This pattern indicates that in the eastern part of the study area the variance in grazing intensity is larger compared to the western part. When plotting the residuals against the fitted values (plot not shown), we see that the large absolute residuals are for low fitted values. Therefore, this may be a consequence of low sample size and the logarithm transformation. I would accept that heterogeneity of variance.

3.2 Influence of grazer density on vegetation traits on a microhabitat scale

dvegc <- dveg[complete.cases(dveg[, c("light.mean", "light.sd", "N.mean", "MV.mean", "horse.corr", "cattle.corr", "year", "cellnr", "heightm", "total_cover")]),]
dvegc$year.z <- as.numeric(scale(dvegc$year))

cor(dvegc[, c("light.mean", "light.sd", "N.mean", "MV.mean", "heightm", "total_cover")])
##              light.mean    light.sd       N.mean      MV.mean     heightm
## light.mean   1.00000000  0.21296431  0.163526150  0.032236377 -0.15073519
## light.sd     0.21296431  1.00000000  0.113114205 -0.040651977  0.01498528
## N.mean       0.16352615  0.11311420  1.000000000 -0.008692668  0.19598873
## MV.mean      0.03223638 -0.04065198 -0.008692668  1.000000000 -0.33577522
## heightm     -0.15073519  0.01498528  0.195988730 -0.335775224  1.00000000
## total_cover -0.06599680  0.12281765  0.020094808  0.029852613  0.33171901
##             total_cover
## light.mean  -0.06599680
## light.sd     0.12281765
## N.mean       0.02009481
## MV.mean      0.02985261
## heightm      0.33171901
## total_cover  1.00000000
min(dvegc$horse.corr[dvegc$horse.corr>0])
## [1] 11.15824
min(dvegc$cattle.corr[dvegc$cattle.corr>0])
## [1] 5.080214
hist(log(dvegc$horse.corr+11))

hist(sqrt(dvegc$horse.corr))

hist(sqrt(dvegc$cattle.corr))

dvegc$horse.corr.sqrt <- sqrt(dvegc$horse.corr)
dvegc$cattle.corr.sqrt <- sqrt(dvegc$cattle.corr)

  
modlightm <- lmer(light.mean ~ horse.corr.sqrt + cattle.corr.sqrt + year.z + (1|cellnr), data=dvegc)
modlightsd <- lmer(light.sd ~ horse.corr.sqrt + cattle.corr.sqrt + year.z + (1|cellnr), data=dvegc)
modN.m <- lmer(N.mean ~ horse.corr.sqrt + cattle.corr.sqrt + year.z + (1|cellnr), data=dvegc)
modMV.m <- lmer(MV.mean ~ horse.corr.sqrt + cattle.corr.sqrt + year.z + (1|cellnr), data=dvegc)
modheightm <- lmer(log(heightm) ~ horse.corr.sqrt + cattle.corr.sqrt + year.z + (1|cellnr), data=dvegc)
modcover <- lmer(total_cover ~ horse.corr.sqrt + cattle.corr.sqrt + year.z + (1|cellnr), data=dvegc)

# plot(modlightm)
# plot(modlightsd)
# plot(modN.m)
# plot(modMV.m)
# plot(modheightm)
# plot(modcover)
# because data is ordered according to cellnr within year, 
# the following plots should show spatial correlation if present
# acf(resid(modlightm))
# acf(resid(modlightsd))
# acf(resid(modN.m))
# acf(resid(modMV.m))
# acf(resid(modheightm))
# acf(resid(modcover))
# -> no spatial correlation

 dvegc$x <- dcells$x[match(dvegc$cellnr, dcells$cellnr)]
 dvegc$y <- dcells$y[match(dvegc$cellnr, dcells$cellnr)]
 
 resdat <- data.frame(residLm=resid(modlightm),
                      residLsd=resid(modlightsd),
                      residNm=resid(modN.m),
                      residMVm=resid(modMV.m),
                      residheight=resid(modheightm),
                      residcover=resid(modcover),
                      x=dvegc$x,
                      y=dvegc$y)
 varLM <- vgram(resdat[,c("x", "y")], resdat$residLm, 
                     type="variogram")
 varLsd <- vgram(resdat[,c("x", "y")], resdat$residLsd, 
                     type="variogram")
 varNm <- vgram(resdat[,c("x", "y")], resdat$residNm, 
                     type="variogram")
 varMVm <- vgram(resdat[,c("x", "y")], resdat$residMVm, 
                     type="variogram")
 varheightm <- vgram(resdat[,c("x", "y")], resdat$residheight, 
                     type="variogram")
 varcover <- vgram(resdat[,c("x", "y")], resdat$residcover, 
                     type="variogram")
 
 
 # par(mfrow=c(3,2))
 # plot(varLM)
 # plot(varLsd)
 # plot(varNm)
 # plot(varMVm)
 # plot(varheightm)
 # plot(varcover)
# -> again, no spatial correlation

nsim <- 5000
bsimLm <- sim(modlightm, n.sim=nsim)
bsimLsd <- sim(modlightsd, n.sim=nsim)
bsimNm <- sim(modN.m, n.sim=nsim)
bsimMVm <- sim(modMV.m, n.sim=nsim)
bsimheightm <- sim(modheightm, n.sim=nsim)
bsimcover <- sim(modcover, n.sim=nsim)

4 Results

4.1 description of the macrohabitat

names(dlandcover)
## [1] "Plot.nr"       "cellnr"        "Cover.tree"    "Cover.sapling"
## [5] "Cover.shrub"   "Cover.meadow"  "Surface.water" "Bareground"
str(dlandcover)
## 'data.frame':    163 obs. of  8 variables:
##  $ Plot.nr      : chr  "A75" "" "A73" "A74" ...
##  $ cellnr       : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ Cover.tree   : int  0 40 0 0 0 30 0 0 0 0 ...
##  $ Cover.sapling: int  3 0 1 2 15 5 1 20 12 1 ...
##  $ Cover.shrub  : int  0 5 0 1 0 0 2 3 3 0 ...
##  $ Cover.meadow : int  97 40 99 97 60 60 97 77 35 69 ...
##  $ Surface.water: int  0 0 0 0 20 5 0 0 50 30 ...
##  $ Bareground   : int  0 15 0 0 5 30 0 0 0 0 ...
landcover_cols <- c("Cover.tree", "Cover.sapling", "Cover.shrub", "Cover.meadow", "Surface.water", "Bareground")
landcover_data <- dlandcover[, landcover_cols]

landcover_data <- as.data.frame(lapply(landcover_data, function(x) as.numeric(as.character(x))))

# convert landcover_data for ggplot2
landcover_long <- melt(landcover_data, variable.name = "LandCoverType", value.name = "Percentage")
## No id variables; using all as measure variables
# mean and sd for each land cover type
landcover_means <- aggregate(Percentage ~ LandCoverType, data = landcover_long, FUN = mean)
landcover_sds <- aggregate(Percentage ~ LandCoverType, data = landcover_long, FUN = sd)

landcover_rename <- c(
  "Cover.tree" = "tree",
  "Cover.sapling" = "sapling",
  "Cover.shrub" = "shrub",
  "Cover.meadow" = "meadow",
  "Surface.water" = "water surface",
  "Bareground" = "bareground"
)

landcover_long$LandCoverType <- landcover_rename[as.character(landcover_long$LandCoverType)]
landcover_means$LandCoverType <- landcover_rename[as.character(landcover_means$LandCoverType)]
landcover_sds$LandCoverType <- landcover_rename[as.character(landcover_sds$LandCoverType)]

landcover_colors <- c(
  "bareground" = "#D6BD8DFF",
  "meadow" = "#F6EDBDFF",
  "sapling" = "#B5B991FF",
  "shrub" = "#778868FF",
  "tree" = "#3D5941FF",
  "water surface" = "#2887A1FF"
)

landcover_long$LandCoverType <- factor(landcover_long$LandCoverType, levels = names(landcover_colors))

# plot 
p <- ggplot(landcover_long, aes(x = LandCoverType, y = Percentage, fill = LandCoverType)) +
  geom_violin(trim = TRUE, scale = "width", adjust = 1.2) +  # adjust the scale and width of the violins
  geom_point(data = landcover_means, aes(x = LandCoverType, y = Percentage), color = "white", size = 3) +
  geom_errorbar(data = landcover_means, 
                    aes(x = LandCoverType, 
                    ymin = pmax(0, Percentage - landcover_sds$Percentage), 
                    ymax = Percentage + landcover_sds$Percentage), 
                    color = "white", width = 0.1) + 
  geom_text(data = landcover_means, aes(x = LandCoverType, y = Percentage, label = round(Percentage, 1)), color = "white", vjust = 1, hjust = -0.4, size = 3.5) +
  scale_fill_manual(values = landcover_colors) +
  labs(title = "Macrohabitat composition on the Rhine Island", x = "Land cover type", y = "Percentage") +
  theme_minimal(base_size = 15) +
  theme(legend.position = "none", panel.background = element_rect(fill = "lightgrey"))

p

4.2 Grazer density in relation to cover variables

4.2.1 Grazers in total (horse and cattle summed)

nsim <- 5000
bsim <- sim(mod, n.sim=nsim)
tab <- t(apply(bsim@fixef, 2, quantile, prob=c(0.5, 0.025, 0.975)))
nfixef <- nrow(tab)


# add random effects
sdcellnr <- apply(bsim@ranef$cellnr[,,1], 1, sd)
tab <- rbind(tab, quantile(sdcellnr, probs=c(0.5, 0.025, 0.975)))
tab <- rbind(tab, quantile(bsim@sigma, probs=c(0.5, 0.025, 0.975)))

rownames(tab)[c(nfixef+1, nfixef+2)] <- c("cellnr_sd", "sigma")

kable(tab, dig=2, caption="Median, 2.5% and 97.5% quantiles of the posterior distributions of the model coefficients for the model with grazer density as response.")
Table 4.1: Median, 2.5% and 97.5% quantiles of the posterior distributions of the model coefficients for the model with grazer density as response.
50% 2.5% 97.5%
(Intercept) -1.81 -2.21 -1.43
Cover.tree 0.00 -0.02 0.02
Cover.sapling 0.00 -0.02 0.03
Cover.shrub -0.09 -0.19 0.02
Surface.water 0.04 0.01 0.08
Bareground 0.02 -0.03 0.07
I(Cover.tree^2) 0.00 0.00 0.00
I(Cover.sapling^2) 0.00 0.00 0.00
I(Cover.shrub^2) 0.01 0.00 0.02
I(Surface.water^2) 0.00 0.00 0.00
I(Bareground^2) 0.00 0.00 0.00
seasonwinter -0.27 -0.53 -0.02
Cover.tree:seasonwinter -0.02 -0.03 0.00
Cover.sapling:seasonwinter 0.02 0.01 0.04
Cover.shrub:seasonwinter 0.18 0.11 0.25
Surface.water:seasonwinter -0.03 -0.05 -0.01
Bareground:seasonwinter 0.00 -0.04 0.03
I(Cover.tree^2):seasonwinter 0.00 0.00 0.00
I(Cover.sapling^2):seasonwinter 0.00 0.00 0.00
I(Cover.shrub^2):seasonwinter -0.02 -0.02 -0.01
I(Surface.water^2):seasonwinter 0.00 0.00 0.00
I(Bareground^2):seasonwinter 0.00 0.00 0.00
cellnr_sd 0.83 0.78 0.88
sigma 0.76 0.72 0.79 
mprop <- apply(dat[,cov.varnames], 2, mean)

par(mfrow=c(3,2), mar=c(4,3,0.5,0.5), mgp=c(2,0.5,0), oma=c(0, 0, 4, 0))
for(cov.var in cov.varnames){
  range.var <- range(dat[,cov.var])
  eval(parse(text=paste0("newdat <- expand.grid(",cov.var,"=seq(range.var[1], range.var[2], length=100),
                         season=levels(factor(dat$season)))")))
  restvar <- cov.varnames[!is.element(cov.varnames, cov.var)]
  for(j in restvar){
    newdat[[j]] <- mprop[j]/(100-mprop[cov.var])*(100-newdat[,cov.var])
  }
  
 #Xmat <- model.matrix(~Cover.tree + Cover.sapling + Cover.shrub+ Surface.water+ Bareground + season + season:(Cover.tree + Cover.sapling + Cover.shrub + Surface.water+ Bareground), data=newdat)
Xmat <- model.matrix(~Cover.tree + Cover.sapling + Cover.shrub+ Surface.water+ Bareground +
              I(Cover.tree^2) + I(Cover.sapling^2) + I(Cover.shrub^2)+ I(Surface.water^2)+ I(Bareground^2) + season + season:(Cover.tree + Cover.sapling + Cover.shrub+ Surface.water+ Bareground + I(Cover.tree^2) + I(Cover.sapling^2) + I(Cover.shrub^2)+ I(Surface.water^2)+ I(Bareground^2)), data=newdat)
  
  
fitmat <- matrix(ncol=nsim, nrow=nrow(newdat))
for(i in 1:nsim) fitmat[,i] <- Xmat%*%bsim@fixef[i,]
newdat$fit <- apply(fitmat, 1, mean)
newdat$lwr <- apply(fitmat, 1, quantile, probs=0.025)
newdat$upr <- apply(fitmat, 1, quantile, probs=0.975)

ylimforlater <- range(dat$grazer.log)
plot(jitter(dat[,cov.var]), dat$grazer.log, xlab="", ylab="log(sum of locations)", col=c("#FFC20A", "#0C7BDC")[as.numeric(dat$season)]) 
title(xlab=paste(cov.var, "[%]"), cex.lab=1.3)

index <- newdat$season=="summer"
lines(newdat[index, cov.var], newdat$fit[index], lwd=2, col="#FFC20A")
lines(newdat[index, cov.var], newdat$lwr[index], lwd=1, lty=3, col="#FFC20A")
lines(newdat[index, cov.var], newdat$upr[index], lwd=1, lty=3, col="#FFC20A")

index <- newdat$season=="winter"
lines(newdat[index, cov.var], newdat$fit[index], lwd=2, col="#0C7BDC")
lines(newdat[index, cov.var], newdat$lwr[index], lwd=1, lty=3, col="#0C7BDC")
lines(newdat[index, cov.var], newdat$upr[index], lwd=1, lty=3, col="#0C7BDC")

mtext("Grazer density", side=3, line=1, outer=TRUE, cex=1.5) 
  
}# close cov.var
Grazer density in relation to the cover variables for summer (orange) and winter (blue). 95% uncertainty intervals are given.

Figure 4.1: Grazer density in relation to the cover variables for summer (orange) and winter (blue). 95% uncertainty intervals are given.

4.2.2 Horse only

nsim <- 5000
bsim <- sim(modh, n.sim=nsim)
tab <- t(apply(bsim@fixef, 2, quantile, prob=c(0.5, 0.025, 0.975)))
nfixef <- nrow(tab)

# add random effects
sdcellnr <- apply(bsim@ranef$cellnr[,,1], 1, sd)
tab <- rbind(tab, quantile(sdcellnr, probs=c(0.5, 0.025, 0.975)))
tab <- rbind(tab, quantile(bsim@sigma, probs=c(0.5, 0.025, 0.975)))

rownames(tab)[c(nfixef+1, nfixef+2)] <- c("cellnr_sd", "sigma")

kable(tab, dig=2, caption="Median, 2.5% and 97.5% quantiles of the posterior distributions of the model coefficients for the model with horse location density as response.")
Table 4.2: Median, 2.5% and 97.5% quantiles of the posterior distributions of the model coefficients for the model with horse location density as response.
50% 2.5% 97.5%
(Intercept) -2.46 -2.84 -2.06
Cover.tree 0.00 -0.02 0.02
Cover.sapling 0.00 -0.02 0.03
Cover.shrub -0.11 -0.22 -0.01
Surface.water 0.02 -0.01 0.05
Bareground 0.02 -0.03 0.07
I(Cover.tree^2) 0.00 0.00 0.00
I(Cover.sapling^2) 0.00 0.00 0.00
I(Cover.shrub^2) 0.01 0.00 0.02
I(Surface.water^2) 0.00 0.00 0.00
I(Bareground^2) 0.00 0.00 0.00
seasonwinter -0.13 -0.37 0.11
Cover.tree:seasonwinter -0.01 -0.02 0.00
Cover.sapling:seasonwinter 0.03 0.01 0.04
Cover.shrub:seasonwinter 0.19 0.13 0.26
Surface.water:seasonwinter 0.00 -0.02 0.02
Bareground:seasonwinter -0.01 -0.04 0.03
I(Cover.tree^2):seasonwinter 0.00 0.00 0.00
I(Cover.sapling^2):seasonwinter 0.00 0.00 0.00
I(Cover.shrub^2):seasonwinter -0.02 -0.02 -0.01
I(Surface.water^2):seasonwinter 0.00 0.00 0.00
I(Bareground^2):seasonwinter 0.00 0.00 0.00
cellnr_sd 0.83 0.79 0.88
sigma 0.69 0.66 0.72 
par(mfrow=c(3,2), mar=c(4,3,0.5,0.5), mgp=c(2,0.5,0), oma=c(0, 0, 4, 0))

for(cov.var in cov.varnames){
  range.var <- range(dat[,cov.var])
  eval(parse(text=paste0("newdat <- expand.grid(",cov.var,"=seq(range.var[1], range.var[2], length=100),
                         season=levels(factor(dat$season)))")))
  restvar <- cov.varnames[!is.element(cov.varnames, cov.var)]
  for(j in restvar){
    newdat[[j]] <- mprop[j]/(100-mprop[cov.var])*(100-newdat[,cov.var])
  }
  
Xmat <- model.matrix(~Cover.tree + Cover.sapling + Cover.shrub+ Surface.water+ Bareground +
              I(Cover.tree^2) + I(Cover.sapling^2) + I(Cover.shrub^2)+ I(Surface.water^2)+ I(Bareground^2) + season + season:(Cover.tree + Cover.sapling + Cover.shrub+ Surface.water+ Bareground + I(Cover.tree^2) + I(Cover.sapling^2) + I(Cover.shrub^2)+ I(Surface.water^2)+ I(Bareground^2)), data=newdat)


fitmat <- matrix(ncol=nsim, nrow=nrow(newdat))
for(i in 1:nsim) fitmat[,i] <- Xmat%*%bsim@fixef[i,]
newdat$fit <- apply(fitmat, 1, mean)
newdat$lwr <- apply(fitmat, 1, quantile, probs=0.025)
newdat$upr <- apply(fitmat, 1, quantile, probs=0.975)

plot(jitter(dat[,cov.var]), dat$horse.log, xlab="", ylab="log(sum of locations)", col=c("#FFC20A", "#0C7BDC")[as.numeric(dat$season)], ylim=ylimforlater)
title(xlab=paste(cov.var, "[%]"), cex.lab=1.3)
index <- newdat$season=="summer"
lines(newdat[index, cov.var], newdat$fit[index], lwd=2, col="#FFC20A")
lines(newdat[index, cov.var], newdat$lwr[index], lwd=1, lty=3, col="#FFC20A")
lines(newdat[index, cov.var], newdat$upr[index], lwd=1, lty=3, col="#FFC20A")

index <- newdat$season=="winter"
lines(newdat[index, cov.var], newdat$fit[index], lwd=2, col="#0C7BDC")
lines(newdat[index, cov.var], newdat$lwr[index], lwd=1, lty=3, col="#0C7BDC")
lines(newdat[index, cov.var], newdat$upr[index], lwd=1, lty=3, col="#0C7BDC")

mtext("Horse density", side=3, line=1, outer=TRUE, cex=1.5)  
}# close cov.var
Horse location density in relation to the cover variables for summer (orange) and winter (blue). 95% uncertainty intervals are given.

Figure 4.2: Horse location density in relation to the cover variables for summer (orange) and winter (blue). 95% uncertainty intervals are given.

4.2.3 Cattle only

nsim <- 5000
bsim <- sim(modc, n.sim=nsim)
tab <- t(apply(bsim@fixef, 2, quantile, prob=c(0.5, 0.025, 0.975)))
nfixef <- nrow(tab)

# add random effects
sdcellnr <- apply(bsim@ranef$cellnr[,,1], 1, sd)
tab <- rbind(tab, quantile(sdcellnr, probs=c(0.5, 0.025, 0.975)))
tab <- rbind(tab, quantile(bsim@sigma, probs=c(0.5, 0.025, 0.975)))

rownames(tab)[c(nfixef+1, nfixef+2)] <- c("cellnr_sd", "sigma")

kable(tab, dig=2, caption="Median, 2.5% and 97.5% quantiles of the posterior distributions of the model coefficients for the model with horse location density as response.")
Table 4.3: Median, 2.5% and 97.5% quantiles of the posterior distributions of the model coefficients for the model with horse location density as response.
50% 2.5% 97.5%
(Intercept) -2.63 -2.98 -2.27
Cover.tree 0.00 -0.02 0.02
Cover.sapling 0.00 -0.02 0.03
Cover.shrub -0.04 -0.13 0.06
Surface.water 0.06 0.03 0.09
Bareground 0.01 -0.03 0.06
I(Cover.tree^2) 0.00 0.00 0.00
I(Cover.sapling^2) 0.00 0.00 0.00
I(Cover.shrub^2) 0.01 0.00 0.01
I(Surface.water^2) 0.00 0.00 0.00
I(Bareground^2) 0.00 0.00 0.00
seasonwinter -0.82 -1.18 -0.47
Cover.tree:seasonwinter -0.01 -0.03 0.01
Cover.sapling:seasonwinter 0.02 0.00 0.05
Cover.shrub:seasonwinter 0.10 0.00 0.19
Surface.water:seasonwinter -0.05 -0.08 -0.02
Bareground:seasonwinter 0.01 -0.04 0.06
I(Cover.tree^2):seasonwinter 0.00 0.00 0.00
I(Cover.sapling^2):seasonwinter 0.00 0.00 0.00
I(Cover.shrub^2):seasonwinter -0.01 -0.02 0.00
I(Surface.water^2):seasonwinter 0.00 0.00 0.00
I(Bareground^2):seasonwinter 0.00 0.00 0.00
cellnr_sd 0.59 0.54 0.64
sigma 1.04 1.00 1.09 
par(mfrow=c(3,2), mar=c(4,3,0.5,0.5), mgp=c(2,0.5,0), oma=c(0, 0, 4, 0))

for(cov.var in cov.varnames){
  range.var <- range(dat[,cov.var])
  eval(parse(text=paste0("newdat <- expand.grid(",cov.var,"=seq(range.var[1], range.var[2], length=100),
                         season=levels(factor(dat$season)))")))
  restvar <- cov.varnames[!is.element(cov.varnames, cov.var)]
  for(j in restvar){
    newdat[[j]] <- mprop[j]/(100-mprop[cov.var])*(100-newdat[,cov.var])
  }
  
Xmat <- model.matrix(~Cover.tree + Cover.sapling + Cover.shrub+ Surface.water+ Bareground + I(Cover.tree^2) + I(Cover.sapling^2) + I(Cover.shrub^2)+ I(Surface.water^2)+ I(Bareground^2) + season + season:(Cover.tree + Cover.sapling + Cover.shrub+ Surface.water+ Bareground + I(Cover.tree^2) + I(Cover.sapling^2) + I(Cover.shrub^2)+ I(Surface.water^2)+ I(Bareground^2)), data=newdat)

fitmat <- matrix(ncol=nsim, nrow=nrow(newdat))
for(i in 1:nsim) fitmat[,i] <- Xmat%*%bsim@fixef[i,]
newdat$fit <- apply(fitmat, 1, mean)
newdat$lwr <- apply(fitmat, 1, quantile, probs=0.025)
newdat$upr <- apply(fitmat, 1, quantile, probs=0.975)


plot(jitter(dat[,cov.var]), dat$cattle.log, xlab="", ylab="log(sum of locations)", col=c("#FFC20A", "#0C7BDC")[as.numeric(dat$season)], ylim=ylimforlater)
title(xlab=paste(cov.var, "[%]"), cex.lab=1.3)

index <- newdat$season=="summer"
lines(newdat[index, cov.var], newdat$fit[index], lwd=2, col="#FFC20A")
lines(newdat[index, cov.var], newdat$lwr[index], lwd=1, lty=3, col="#FFC20A")
lines(newdat[index, cov.var], newdat$upr[index], lwd=1, lty=3, col="#FFC20A")

index <- newdat$season=="winter"
lines(newdat[index, cov.var], newdat$fit[index], lwd=2, col="#0C7BDC")
lines(newdat[index, cov.var], newdat$lwr[index], lwd=1, lty=3, col="#0C7BDC")
lines(newdat[index, cov.var], newdat$upr[index], lwd=1, lty=3, col="#0C7BDC")

mtext("Cattle density", side=3, line=1, outer=TRUE, cex=1.5)

}# close cov.var
Cattle location density in relation to the cover variables for summer (orange) and winter (blue). 95% uncertainty intervals are given.

Figure 4.3: Cattle location density in relation to the cover variables for summer (orange) and winter (blue). 95% uncertainty intervals are given.

4.3 Relationship between grazer density and vegetation traits on the microhabitat scale

Diversity of vegetation on the study site in the different years. The diversity indices are measured taking the number of occupied cells per species as if they were number of individuals.

dveglong <- data.frame(species=c(dveg$species1, dveg$species2, dveg$species3),
                       year=rep(dveg$year, 3),
                       season=rep(dveg$season,3),
                       plot_nr=rep(dveg$plot_nr,3))
ddiv <- aggregate(plot_nr~species + year+ season, data=dveglong, FUN=length)

ddiversity <- data.frame(year=2018:2021,
                         nrspecies = as.numeric(table(ddiv$year)),
                         totalcells=tapply(ddiv$plot_nr, ddiv$year, FUN=sum))
ddiv$totalcells <- ddiversity$totalcells[match(ddiv$year, ddiversity$year)]
ddiv$p <- ddiv$plot_nr/ddiv$totalcells
shannon <- function(p) -sum(p*log(p)) # log is natural logarithm
ddiversity$Shannon <- tapply(ddiv$p, ddiv$year, FUN=shannon)
ddiversity$Pielou_evenness <- ddiversity$Shannon/log(ddiversity$nrspecies)

kable(ddiversity[,-3], dig=3, caption="Species richness (nrspecies), Shannon diversity index and Pielou evenness for each year")
Table 4.4: Species richness (nrspecies), Shannon diversity index and Pielou evenness for each year
year nrspecies Shannon Pielou_evenness
2018 2018 51 3.352321 0.8526119
2019 2019 46 3.263625 0.8524238
2020 2020 42 3.099439 0.8292438
2021 2021 43 3.138654 0.8344821
Table 4.5: Estimated mean (Intercept) and partial coefficients of grazing (horse and cattle) and year. Year was standardised before model fitting.
Interc. Int.lwr Int.upr horse horse.lwr horse.upr cattle cattle.lwr cattle.upr year.z year.lwr year.upr
Light.mean 3.38 3.30 3.46 0.00 0.00 0.01 0.00 0.00 0.01 -0.03 -0.07 0.01
Light.SD 0.40 0.33 0.47 0.00 -0.01 0.01 0.01 0.00 0.02 -0.02 -0.06 0.01
N.mean 3.02 2.90 3.15 0.00 -0.02 0.01 0.01 -0.01 0.02 -0.07 -0.13 -0.01
MV.mean 2.89 2.73 3.04 0.01 0.00 0.03 -0.02 -0.04 0.00 0.23 0.15 0.31
log(height.m) 2.45 2.28 2.62 -0.01 -0.02 0.01 -0.01 -0.03 0.01 -0.25 -0.34 -0.17
total cover 70.65 64.63 76.72 0.67 0.15 1.22 -0.30 -0.85 0.24 -0.06 -2.55 2.44 
newdathorse <- data.frame(horse.corr=seq(0, 510, length=100),
                          cattle.corr=mean(dvegc$cattle.corr),
                          year.z=0)
newdathorse$horse.corr.sqrt <- sqrt(newdathorse$horse.corr)
newdathorse$cattle.corr.sqrt <- sqrt(newdathorse$cattle.corr)

newdatcattle <- data.frame(cattle.corr=seq(0, 800, length=100),
                          horse.corr=mean(dvegc$horse.corr),
                          year.z=0)
newdatcattle$horse.corr.sqrt <- sqrt(newdatcattle$horse.corr)
newdatcattle$cattle.corr.sqrt <- sqrt(newdatcattle$cattle.corr)

newdatyear <- data.frame(cattle.corr=mean(dvegc$cattle.corr),
                          horse.corr=mean(dvegc$horse.corr),
                          year=2018:2021)
newdatyear$horse.corr.sqrt <- sqrt(newdatyear$horse.corr)
newdatyear$cattle.corr.sqrt <- sqrt(newdatyear$cattle.corr)
newdatyear$year.z <- (newdatyear$year-mean(dvegc$year))/sd(dvegc$year)

Xmathorse <- model.matrix(~horse.corr.sqrt + cattle.corr.sqrt + year.z, data=newdathorse)
Xmatcattle <- model.matrix(~horse.corr.sqrt + cattle.corr.sqrt + year.z, data=newdatcattle)
Xmatyear <- model.matrix(~horse.corr.sqrt + cattle.corr.sqrt + year.z, data=newdatyear)

# Light mean
newdathorse$Lm <- predict(modlightm, newdata=newdathorse, re.form=NA)
fitmat <- matrix(ncol=nsim, nrow=nrow(newdathorse))
for(i in 1:nsim) fitmat[,i] <- Xmathorse%*%bsimLm@fixef[i,]
newdathorse$Lm.lwr <- apply(fitmat, 1, quantile, prob=0.025)
newdathorse$Lm.upr <- apply(fitmat, 1, quantile, prob=0.975)

newdatcattle$Lm <- predict(modlightm, newdata=newdatcattle, re.form=NA)
fitmat <- matrix(ncol=nsim, nrow=nrow(newdatcattle))
for(i in 1:nsim) fitmat[,i] <- Xmatcattle%*%bsimLm@fixef[i,]
newdatcattle$Lm.lwr <- apply(fitmat, 1, quantile, prob=0.025)
newdatcattle$Lm.upr <- apply(fitmat, 1, quantile, prob=0.975)

newdatyear$Lm <- predict(modlightm, newdata=newdatyear, re.form=NA)
fitmat <- matrix(ncol=nsim, nrow=nrow(newdatyear))
for(i in 1:nsim) fitmat[,i] <- Xmatyear%*%bsimLm@fixef[i,]
newdatyear$Lm.lwr <- apply(fitmat, 1, quantile, prob=0.025)
newdatyear$Lm.upr <- apply(fitmat, 1, quantile, prob=0.975)

# Light SD
newdathorse$Lsd <- predict(modlightsd, newdata=newdathorse, re.form=NA)
fitmat <- matrix(ncol=nsim, nrow=nrow(newdathorse))
for(i in 1:nsim) fitmat[,i] <- Xmathorse%*%bsimLsd@fixef[i,]
newdathorse$Lsd.lwr <- apply(fitmat, 1, quantile, prob=0.025)
newdathorse$Lsd.upr <- apply(fitmat, 1, quantile, prob=0.975)

newdatcattle$Lsd <- predict(modlightsd, newdata=newdatcattle, re.form=NA)
fitmat <- matrix(ncol=nsim, nrow=nrow(newdatcattle))
for(i in 1:nsim) fitmat[,i] <- Xmatcattle%*%bsimLsd@fixef[i,]
newdatcattle$Lsd.lwr <- apply(fitmat, 1, quantile, prob=0.025)
newdatcattle$Lsd.upr <- apply(fitmat, 1, quantile, prob=0.975)

newdatyear$Lsd <- predict(modlightsd, newdata=newdatyear, re.form=NA)
fitmat <- matrix(ncol=nsim, nrow=nrow(newdatyear))
for(i in 1:nsim) fitmat[,i] <- Xmatyear%*%bsimLsd@fixef[i,]
newdatyear$Lsd.lwr <- apply(fitmat, 1, quantile, prob=0.025)
newdatyear$Lsd.upr <- apply(fitmat, 1, quantile, prob=0.975)

# N mean
newdathorse$Nm <- predict(modN.m, newdata=newdathorse, re.form=NA)
fitmat <- matrix(ncol=nsim, nrow=nrow(newdathorse))
for(i in 1:nsim) fitmat[,i] <- Xmathorse%*%bsimNm@fixef[i,]
newdathorse$Nm.lwr <- apply(fitmat, 1, quantile, prob=0.025)
newdathorse$Nm.upr <- apply(fitmat, 1, quantile, prob=0.975)

newdatcattle$Nm <- predict(modN.m, newdata=newdatcattle, re.form=NA)
fitmat <- matrix(ncol=nsim, nrow=nrow(newdatcattle))
for(i in 1:nsim) fitmat[,i] <- Xmatcattle%*%bsimNm@fixef[i,]
newdatcattle$Nm.lwr <- apply(fitmat, 1, quantile, prob=0.025)
newdatcattle$Nm.upr <- apply(fitmat, 1, quantile, prob=0.975)

newdatyear$Nm <- predict(modN.m, newdata=newdatyear, re.form=NA)
fitmat <- matrix(ncol=nsim, nrow=nrow(newdatyear))
for(i in 1:nsim) fitmat[,i] <- Xmatyear%*%bsimNm@fixef[i,]
newdatyear$Nm.lwr <- apply(fitmat, 1, quantile, prob=0.025)
newdatyear$Nm.upr <- apply(fitmat, 1, quantile, prob=0.975)

#MV mean
newdathorse$MVm <- predict(modMV.m, newdata=newdathorse, re.form=NA)
fitmat <- matrix(ncol=nsim, nrow=nrow(newdathorse))
for(i in 1:nsim) fitmat[,i] <- Xmathorse%*%bsimMVm@fixef[i,]
newdathorse$MVm.lwr <- apply(fitmat, 1, quantile, prob=0.025)
newdathorse$MVm.upr <- apply(fitmat, 1, quantile, prob=0.975)

newdatcattle$MVm <- predict(modMV.m, newdata=newdatcattle, re.form=NA)
fitmat <- matrix(ncol=nsim, nrow=nrow(newdatcattle))
for(i in 1:nsim) fitmat[,i] <- Xmatcattle%*%bsimMVm@fixef[i,]
newdatcattle$MVm.lwr <- apply(fitmat, 1, quantile, prob=0.025)
newdatcattle$MVm.upr <- apply(fitmat, 1, quantile, prob=0.975)

newdatyear$MVm <- predict(modMV.m, newdata=newdatyear, re.form=NA)
fitmat <- matrix(ncol=nsim, nrow=nrow(newdatyear))
for(i in 1:nsim) fitmat[,i] <- Xmatyear%*%bsimMVm@fixef[i,]
newdatyear$MVm.lwr <- apply(fitmat, 1, quantile, prob=0.025)
newdatyear$MVm.upr <- apply(fitmat, 1, quantile, prob=0.975)


#height mean
newdathorse$heightm <- predict(modheightm, newdata=newdathorse, re.form=NA)
fitmat <- matrix(ncol=nsim, nrow=nrow(newdathorse))
for(i in 1:nsim) fitmat[,i] <- Xmathorse%*%bsimheightm@fixef[i,]
newdathorse$heightm.lwr <- apply(fitmat, 1, quantile, prob=0.025)
newdathorse$heightm.upr <- apply(fitmat, 1, quantile, prob=0.975)

newdatcattle$heightm <- predict(modheightm, newdata=newdatcattle, re.form=NA)
fitmat <- matrix(ncol=nsim, nrow=nrow(newdatcattle))
for(i in 1:nsim) fitmat[,i] <- Xmatcattle%*%bsimheightm@fixef[i,]
newdatcattle$heightm.lwr <- apply(fitmat, 1, quantile, prob=0.025)
newdatcattle$heightm.upr <- apply(fitmat, 1, quantile, prob=0.975)

newdatyear$heightm <- predict(modheightm, newdata=newdatyear, re.form=NA)
fitmat <- matrix(ncol=nsim, nrow=nrow(newdatyear))
for(i in 1:nsim) fitmat[,i] <- Xmatyear%*%bsimheightm@fixef[i,]
newdatyear$heightm.lwr <- apply(fitmat, 1, quantile, prob=0.025)
newdatyear$heightm.upr <- apply(fitmat, 1, quantile, prob=0.975)

#cover mean (%)
newdathorse$cover <- predict(modcover, newdata=newdathorse, re.form=NA)
fitmat <- matrix(ncol=nsim, nrow=nrow(newdathorse))
for(i in 1:nsim) fitmat[,i] <- Xmathorse%*%bsimcover@fixef[i,]
newdathorse$cover.lwr <- apply(fitmat, 1, quantile, prob=0.025)
newdathorse$cover.upr <- apply(fitmat, 1, quantile, prob=0.975)

newdatcattle$cover <- predict(modcover, newdata=newdatcattle, re.form=NA)
fitmat <- matrix(ncol=nsim, nrow=nrow(newdatcattle))
for(i in 1:nsim) fitmat[,i] <- Xmatcattle%*%bsimcover@fixef[i,]
newdatcattle$cover.lwr <- apply(fitmat, 1, quantile, prob=0.025)
newdatcattle$cover.upr <- apply(fitmat, 1, quantile, prob=0.975)

newdatyear$cover <- predict(modcover, newdata=newdatyear, re.form=NA)
fitmat <- matrix(ncol=nsim, nrow=nrow(newdatyear))
for(i in 1:nsim) fitmat[,i] <- Xmatyear%*%bsimcover@fixef[i,]
newdatyear$cover.lwr <- apply(fitmat, 1, quantile, prob=0.025)
newdatyear$cover.upr <- apply(fitmat, 1, quantile, prob=0.975)

#light mean
par(mfrow=c(6,3), mar=c(2,2,0.5,0.5), oma=c(0,2,2,0))
plot(dvegc$horse.corr, dvegc$light.mean, pch=16, col="#2C6B4F66", cex=0.5, main=NA, ylab="Light mean", xpd=NA, xlab=NA)
mtext("Horse density", side=3, line=0.5)
lines(newdathorse$horse.corr, newdathorse$Lm, lwd=2)
lines(newdathorse$horse.corr, newdathorse$Lm.lwr, lwd=1,lty=3)
lines(newdathorse$horse.corr, newdathorse$Lm.upr, lwd=1,lty=3)

plot(dvegc$cattle.corr, dvegc$light.mean, pch=16, col="#2C6B4F66", cex=0.5, main=NA)
mtext("Cattle density", side=3, line=0.5)
lines(newdatcattle$cattle.corr, newdatcattle$Lm, lwd=2)
lines(newdatcattle$cattle.corr, newdatcattle$Lm.lwr, lwd=1,lty=3)
lines(newdatcattle$cattle.corr, newdatcattle$Lm.upr, lwd=1,lty=3)

plot(jitter(dvegc$year), dvegc$light.mean, pch=16, col="#2C6B4F66", cex=0.5, main=NA)
mtext("Year", side=3, line=0.5)
lines(newdatyear$year, newdatyear$Lm, lwd=2)
lines(newdatyear$year, newdatyear$Lm.lwr, lwd=1,lty=3)
lines(newdatyear$year, newdatyear$Lm.upr, lwd=1,lty=3)

# light sd

plot(dvegc$horse.corr, dvegc$light.sd, pch=16, col="#2C6B4F66", cex=0.5, ylab="Light SD", xpd=NA, xlab="")
lines(newdathorse$horse.corr, newdathorse$Lsd, lwd=2)
lines(newdathorse$horse.corr, newdathorse$Lsd.lwr, lwd=1,lty=3)
lines(newdathorse$horse.corr, newdathorse$Lsd.upr, lwd=1,lty=3)

plot(dvegc$cattle.corr, dvegc$light.sd, pch=16, col="#2C6B4F66", cex=0.5)
lines(newdatcattle$cattle.corr, newdatcattle$Lsd, lwd=2)
lines(newdatcattle$cattle.corr, newdatcattle$Lsd.lwr, lwd=1,lty=3)
lines(newdatcattle$cattle.corr, newdatcattle$Lsd.upr, lwd=1,lty=3)

plot(jitter(dvegc$year), dvegc$light.sd, pch=16, col="#2C6B4F66", cex=0.5)
lines(newdatyear$year, newdatyear$Lsd, lwd=2)
lines(newdatyear$year, newdatyear$Lsd.lwr, lwd=1,lty=3)
lines(newdatyear$year, newdatyear$Lsd.upr, lwd=1,lty=3)

# N mean
plot(dvegc$horse.corr, dvegc$N.mean, pch=16, col="#2C6B4F66", cex=0.5, ylab="N mean", xpd=NA, xlab=NA)
lines(newdathorse$horse.corr, newdathorse$Nm, lwd=2)
lines(newdathorse$horse.corr, newdathorse$Nm.lwr, lwd=1,lty=3)
lines(newdathorse$horse.corr, newdathorse$Nm.upr, lwd=1,lty=3)

plot(dvegc$cattle.corr, dvegc$N.mean, pch=16, col="#2C6B4F66", cex=0.5)
lines(newdatcattle$cattle.corr, newdatcattle$Nm, lwd=2)
lines(newdatcattle$cattle.corr, newdatcattle$Nm.lwr, lwd=1,lty=3)
lines(newdatcattle$cattle.corr, newdatcattle$Nm.upr, lwd=1,lty=3)

plot(jitter(dvegc$year), dvegc$N.mean, pch=16, col="#2C6B4F66", cex=0.5)
lines(newdatyear$year, newdatyear$Nm, lwd=2)
lines(newdatyear$year, newdatyear$Nm.lwr, lwd=1,lty=3)
lines(newdatyear$year, newdatyear$Nm.upr, lwd=1,lty=3)

# MV mean
plot(dvegc$horse.corr, dvegc$MV.mean, pch=16, col="#2C6B4F66", cex=0.5, ylab="MV mean", xpd=NA, xlab="")
lines(newdathorse$horse.corr, newdathorse$MVm, lwd=2)
lines(newdathorse$horse.corr, newdathorse$MVm.lwr, lwd=1,lty=3)
lines(newdathorse$horse.corr, newdathorse$MVm.upr, lwd=1,lty=3)

plot(dvegc$cattle.corr, dvegc$MV.mean, pch=16, col="#2C6B4F66", cex=0.5)
lines(newdatcattle$cattle.corr, newdatcattle$MVm, lwd=2)
lines(newdatcattle$cattle.corr, newdatcattle$MVm.lwr, lwd=1,lty=3)
lines(newdatcattle$cattle.corr, newdatcattle$MVm.upr, lwd=1,lty=3)

plot(jitter(dvegc$year), dvegc$MV.mean, pch=16, col="#2C6B4F66", cex=0.5)
lines(newdatyear$year, newdatyear$MVm, lwd=2)
lines(newdatyear$year, newdatyear$MVm.lwr, lwd=1,lty=3)
lines(newdatyear$year, newdatyear$MVm.upr, lwd=1,lty=3)

# height m
plot(dvegc$horse.corr, dvegc$heightm, pch=16, col="#2C6B4F66", cex=0.5, ylab="Height mean", xpd=NA, xlab="")
lines(newdathorse$horse.corr, exp(newdathorse$heightm), lwd=2)
lines(newdathorse$horse.corr, exp(newdathorse$heightm.lwr), lwd=1,lty=3)
lines(newdathorse$horse.corr, exp(newdathorse$heightm.upr), lwd=1,lty=3)

plot(dvegc$cattle.corr, dvegc$heightm, pch=16, col="#2C6B4F66", cex=0.5)
lines(newdatcattle$cattle.corr, exp(newdatcattle$heightm), lwd=2)
lines(newdatcattle$cattle.corr, exp(newdatcattle$heightm.lwr), lwd=1,lty=3)
lines(newdatcattle$cattle.corr, exp(newdatcattle$heightm.upr), lwd=1,lty=3)

plot(jitter(dvegc$year), dvegc$heightm, pch=16, col="#2C6B4F66", cex=0.5)
lines(newdatyear$year, exp(newdatyear$heightm), lwd=2)
lines(newdatyear$year, exp(newdatyear$heightm.lwr), lwd=1,lty=3)
lines(newdatyear$year, exp(newdatyear$heightm.upr), lwd=1,lty=3)

# cover
plot(dvegc$horse.corr, dvegc$total_cover, pch=16, col="#2C6B4F66", cex=0.5, ylab="Total cover", xpd=NA, xlab="")
lines(newdathorse$horse.corr, newdathorse$cover, lwd=2)
lines(newdathorse$horse.corr, newdathorse$cover.lwr, lwd=1,lty=3)
lines(newdathorse$horse.corr, newdathorse$cover.upr, lwd=1,lty=3)

plot(dvegc$cattle.corr, dvegc$total_cover, pch=16, col="#2C6B4F66", cex=0.5)
lines(newdatcattle$cattle.corr, newdatcattle$cover, lwd=2)
lines(newdatcattle$cattle.corr, newdatcattle$cover.lwr, lwd=1,lty=3)
lines(newdatcattle$cattle.corr, newdatcattle$cover.upr, lwd=1,lty=3)

plot(jitter(dvegc$year), dvegc$total_cover, pch=16, col="#2C6B4F66", cex=0.5)
lines(newdatyear$year, newdatyear$cover, lwd=2)
lines(newdatyear$year, newdatyear$cover.lwr, lwd=1,lty=3)
lines(newdatyear$year, newdatyear$cover.upr, lwd=1,lty=3)
Relationship between horse and cattle densities, and plant functional and structural traits

Figure 4.4: Relationship between horse and cattle densities, and plant functional and structural traits